function dataset = generate_counterfactual(Prob, Ss, rs, Hs, tol, max_It, rs_series, type);
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha...
        xi xi1 ub lb pi_e H w theta Int_pi_w...
        sigmae elasts Pr log_e_mu log_e_sigma skill_shock  phi sigma gamma rts fp r_bar r_base w_mult shares...
        pe thetas tmp_mult adj w_bar ss_calc unemp h_resid shock_multiplier myslope eta mu_a thetas
    
    eta = 1;
   
    % Empty matrix for data
    dataset = cell(5, 4);
    ln_e_hat = -1/2*log_e_sigma^2; % mean of the non-normalized idiosyncratic wage shock
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial Choices
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
    
    educ_mat_ss=educ_mat;
    myDist0 = P_S .* Prob;
    myDist0  = sum(bsxfun(@times, myDist0 , reshape(pi_theta, [1,1,1,2])),4);
    N_size = squeeze(sum(myDist0,1:2));
    avg_ability = squeeze(sum(a'.* myDist0,1:2)./sum( myDist0 ,1:2))';
    avg_logability = squeeze(sum(log(a)'.* myDist0,1:2)./sum( myDist0 ,1:2))';

    N_n = length(Ss);
    
    denom = sum(myDist0 .* wprimeS, 1:3);
    bet1 = bet1 / (denom)^xi1;
    Ss = Ss .* (denom)^xi1;
    
    Hs =  inv_hs(shares , rs); 

    w_max0 = max(w);
    w_min0 = min(w);

    Probm = update_dist_dyn(Prob, Pr, Ss, rs, 1);
    F_es = skill_shock .* F_es;
    b0 =  b ;
   
    Prob = Probm;

    myDist0 = P_S .* Prob;
    myDist0  = sum(bsxfun(@times, myDist0 , reshape(pi_theta, [1,1,1,2])),4);

    for t = 1:17
        
        if t == 1
           eta = skill_shock;
        end
        pi_w = sum(myDist0,2:3);
        w_bar = dot(w, pi_w);

        w0 = w;
        
        if type == 1

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Random Allocation
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            r_0 = rs(end);

            for liln = 1:12220
              f_0 = (1.0893)^t - sum((Hs.* r_0).^elasts);
              f_p_0 = -dot(elasts, (Hs .* r_0 ).^(elasts-1) .* Hs );
              r_1 = r_0 - f_0 / f_p_0 *0.1;
              eps = abs(r_1 - r_0);
              r_0 = r_1;
              if eps < 1e-6
                  break
              elseif liln == 12220
                  error("Did not clear housing market")
              end
            end
            rs = ones(length(N_size),1)*r_0;
            % Random allocation proportional to N size
            P_S = ones(size(P_S)) .* reshape(N_size, [1 1 3 1]);
            
            for iter = 1:540
                [~, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
                P_S2 = sum(bsxfun(@times, P_S, reshape(pi_theta, [1,1,1,2])),4);
                popMass = bsxfun(@times, P_S2, Prob);
                wMass = squeeze(sum(sum(bsxfun(@times, wprimeS, popMass), 1:2), 4));
                popMass = squeeze(sum(sum(popMass, 1:2),4));
                
                Ssi = wMass ./ popMass;
                Ssi = Ssi(:).^xi1; 
                d_S = abs(Ssi - Ss);
                if sum(d_S) < 2 * tol
                    break
                end
                Ss = Ssi;
            end       
            % Future dist
            Probi = update_dist_int_edu_rm(Prob, Pr, Ss, rs, N_size, 0);
            % Current dist
            myDist1 = P_S .* Prob;
            myDist1  = sum(bsxfun(@times, myDist1 , reshape(pi_theta, [1,1,1,2])),4);
        elseif type == 2
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % SS allocations - fixed price
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           
            for iter = 1:540
                [~, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
                wMass = squeeze(sum(sum(bsxfun(@times, wprimeS, myDist0), 1:2), 4));

                popMass = squeeze(sum(sum(myDist0, 1:2),4));
                Ssi = wMass ./ popMass;
                Ssi = Ssi(:).^xi1;
                d_S = abs(Ssi - Ss);
                if sum(d_S) < 2 * tol / 10
                    break
                end
                Ss = Ssi;
            end
            
            myDist1 = myDist0;
            myDist0 = update_dist_int_edu_ss(myDist0, Pr, Ss, rs);
            Probi = squeeze(sum(myDist0,3));
        
        elseif type == 3
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % SS allocations - prices adjust
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            rs = ((N_size*(1.0893^t)).^(1./elasts) )./Hs ;
            
            for iter = 1:40
                [~, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
                wMass = squeeze(sum(sum(bsxfun(@times, wprimeS, myDist0), 1:2), 4));

                popMass = squeeze(sum(sum(myDist0, 1:2),4));
                Ssi = wMass ./ popMass;
                Ssi = Ssi(:).^xi1;
                d_S = abs(Ssi - Ss);
                if sum(d_S) < 2 * tol
                    break
                end
                Ss = Ssi;
            end
            
            myDist1 = myDist0;
            myDist0 = update_dist_int_edu_ss(myDist0, Pr, Ss, rs);
            Probi = squeeze(sum(myDist0,3));

        
        elseif type==4
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Fixed location and education
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
     
             for iter = 1:540
                [~, F_es, wprimeS, educ_mat_ss] = choiceProb_no_ed(Ss, rs,educ_mat_ss);
                wMass = squeeze(sum(sum(bsxfun(@times, wprimeS, myDist0), 1:2), 4));

                popMass = squeeze(sum(sum(myDist0, 1:2),4));
                Ssi = wMass ./ popMass;
                Ssi = Ssi(:).^xi1;
                d_S = abs(Ssi - Ss);
                if sum(d_S) < 2 * tol / 10
                    break
                end
                Ss = Ssi;
            end
            
            myDist1 = myDist0;
            myDist0 = update_dist_noedu(myDist0, Pr, Ss, rs,educ_mat_ss);
            Probi = squeeze(sum(myDist0,3));
        end
        Prob_kids = Probi;
        pi_w = squeeze(sum(myDist1, 2:3));
        cdf_w = cumsum(pi_w,1);
        avg_ed = squeeze(sum(educ_mat.* myDist1,1:2)./sum( myDist1 ,1:2))';
        avg_ability = squeeze(sum(a' .* myDist1,1:2)./sum( myDist1 ,1:2))';
        avg_logability = squeeze(sum(log(a)' .* myDist1,1:2)./sum( myDist1 ,1:2))';
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Gini Index
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        my_prod = pi_w .* pi_w' .* abs(w - w');
        gini = sum(my_prod(:))/(2*dot(pi_w, w));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Dissimilarity Index
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        rich_cutoff = 0.8;
        poor = rich_cutoff(1); 
        rich = 1 - rich_cutoff(1); 
        agg_cdf = cumsum(myDist1, 1);
       
        find_rich_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-rich_cutoff;
        w_rich = fzero(find_rich_w,(w_max0+w_min0)/2);
        
        share_h = squeeze(sum((w0 < w_rich) .* myDist1,1:2))/ sum(squeeze(sum((w0 < w_rich) .* myDist1,1:2)));
        share_l = (squeeze(sum((w0 >= w_rich) .* myDist1,1:2)))/(sum(squeeze(sum((w0 >= w_rich) .* myDist1,1:2))));
        
        dissim = sum(abs(share_h - share_l))*0.5; 
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Rank-Rank Correlation
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Upward Mobility
        % find 25th percentile (cdf_w was already defined above)
        N_bp1 = 125;
        q = linspace(0,1,N_bp1);
        w_q = zeros(1,N_bp1);
        w_q(1) = w_min0;
        w_q(N_bp1)= w_max0+1e-5;
        for i = 1:(N_bp1-2)
          find_q_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-q(i+1);
          w_q(i+1) = fzero(find_q_w,(w_max0+w_min0)/2);
        end
        
        w_qk = zeros(1,N_bp1);
        w_qk(1) = w_min;
        w_qk(N_bp1)= w_max+1e-5;
        cdf_w_kids = cumsum(sum(Prob_kids,2));
        for i = 1:(N_bp1-2)
          find_q_w = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-q(i+1);
          w_qk(i+1) = fzero(find_q_w,(w_max+w_min)/2);
        end
        
        parents_rank = zeros(N_w,N_a);
        for m = 1:(N_bp1-1)
            for i = 1:N_w
                if w_q(m)<=w0(i) && w_q(m+1)>w0(i)
                    parents_rank(i,:) = m;
                end
            end
        end

        wprime = unemp + (b0 + F_es) .* f_w(w0);
        wprime = repmat(wprime, [1 1 1 N_e 2]);
        tmp = repmat(e, [1 N_w N_a N_n 2]);
        tmp = permute(tmp, [2 3 4 1 5]);
       
        wprime = wprime .* tmp;
        
        [~, kids_rank] = histc(wprime,w_qk);
        pi_joint_w_a_e = repmat(myDist1,1,1,1,N_e) .* permute(repmat(pi_e, [1 N_w N_a N_n]), [2 3 4 1]);
        mu_prank = sum(sum(sum(sum(parents_rank.*pi_joint_w_a_e))));
        mu_krank = sum(sum(sum(sum(kids_rank.*pi_joint_w_a_e))));
        mu_krank = pi_theta*mu_krank(:);
        prank_std = sum(sum(sum(sum((parents_rank-mu_prank).^2.*pi_joint_w_a_e))));
        prank_std = sqrt(prank_std(:));
        krank_std = sum(sum(sum(sum((kids_rank-mu_krank).^2.*pi_joint_w_a_e))));
        krank_std = sqrt(pi_theta*krank_std(:));
        rank_cov = sum(sum(sum(sum((parents_rank - mu_prank).*(kids_rank-mu_krank).*pi_joint_w_a_e))));
        rank_cov = pi_theta*rank_cov(:);

        rank_corr = rank_cov/(prank_std*krank_std);
       

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Results
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        dataset{t,1} = dissim;
        dataset{t,2} = gini;
        dataset{t,3} = Ss';        
        dataset{t,4} = rs'; 
        dataset{t,5} = avg_ability; 
        dataset{t,6} = rank_corr; 
        dataset{t,7} = avg_logability;
        
        if ( type == 1) | (type == 4) 
            Prob = Probi;
          
        end
       
    end

end